• O R é um ambiente de computação estatística, com mais de 13 mil pacotes publicados.
• Cada um destes pacotes tem um fim específico. Neste curso, utilizaremos principalmente as bibliotecas ‘forecast’ e o fpp2. Ambos os pacotes são utlizados no livro Forecating: Principles and Practice, bem como os exemplos destes slides também são baseados no trabalho do autor deste livro.

• Você trabalha em uma empresa e a mesma não utiliza muitas técnicas para a previsão de demanda do produto “x”.

• Após estudar previsão de demanda, você percebeu que as vendas deste produto possuem uma sazonalidade anual (o que é óbvio), porém a empresa não utiliza técnica de previsão de demanda.

• Se a empresa possuísse uma previsão melhor, seria possível atuar de modo mais estratégico sobre o planejamento da produção (talvez antecipar a produção em meses de pico ao invés de realizar horas extra).

• Você já pesquisou alguns softwares para previsão de demanda, porém ou são complicados de usar, ou possuem um preço proibitivo.

1. Faça um gráfico dos dados;
2. Selecione uma função para a previsão;
3. Estime os parâmetros da função;
4. Avalie a Qualidade do modelo de previsão;
5. Selecione e implemente o melhor modelo encontrado.

• Instale E carregue as bibliotecas que iremos utilizar nesta aula rodando estes comandos:

bibliotecas = c("forecast", "fpp2", "readxl")

install.packages(bibliotecas)

# Caso encontre algum erro, rode o install packages separadamente para cada biblioteca:
install.packages("forecast")
install.packages("fpp2")
install.packages("readxl")

library(forecast)
library(fpp2)
library(readxl)

• Observe as série gold:

# A série temporal "gold" foi carregada pela biblioteca forecast.
autoplot(gold)

• Tendência: Os dados possuem uma tendência geral de aumento ou queda (exemplo: Crescimento da população).
• Sazonalidade: Os dados se comportam de modo similar, obedecendo padrões com durações fixas (exemplo: venda de ovos de páscoa).
• Ciclicidade: A série possui um comportamento variando em ciclos, porém sem duração fixa (exemplo: ciclos econômicos.).

# Produção de Lâ na Austrália
autoplot(woolyrnq) 

# Produção de Gás na Austrália
autoplot(gas) 

# A função frequency determina a frequência da série. 
# Para dados sazonais, irá definir o período dominante da sazonalidade, 
# e para dados em ciclos, a duração média dos ciclos.
frequency(gas) 

## [1] 12

Observando a série temporal:

# Produção de Gás na Austrália
autoplot(a10) + ylab("Demanda Anti-Diabeticos")

Observando um gráfico sazonal

ggseasonplot(a10) 

Observando um gráfico sazonal “polar”:

ggseasonplot(a10, polar = T)

Observando a Venda de Cerveja:

beer <- window(ausbeer, start = 1992)
autoplot(beer) + ylab("Venda de Cerveja")

Observando a Venda de Cerveja e sua Sazonalidade:

ggseasonplot(beer) 

Mais uma forma de visualizar a demanda sazonal

ggsubseriesplot(beer)

• Para manipular séries temporais no R, é interessante criar objetos do tipo ‘ts’.
• Podemos criar objetos a partir de arquivos do excel, csv, bases de dados, ou mesmo APIs públicas.

• Podemos utilizar a biblioteca ‘readxl’ para ler arquivos do excel:

library(readxl)
dados_excel = readxl::read_xlsx(path = "dados/VendasCarros.xlsx", sheet = "Dados")
head(dados_excel)

## # A tibble: 6 x 4
##     Ano   Mês VendasCarros VendasCarrosAux
##   <dbl> <dbl>        <dbl>           <dbl>
## 1  2018     1         2493            2486
## 2  2018     2         2875            3003
## 3  2018     3         3504            3454
## 4  2018     4         3786            3665
## 5  2018     5         4094            4353
## 6  2018     6         4994            4842

• Para trabalhar com séries temporais, iremos transformar esta tabela:

ts_vendas = ts(data = dados_excel$VendasCarros, start = c(2018,1),frequency = 12)
ts_vendas

##       Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec
## 2018 2493 2875 3504 3786 4094 4994 4910 5575 5839 6202 6833 7382
## 2019 2470 3059 3537 4003 4110 4516 5364 5854 5752 6089 7025 7801
## 2020 2473 2884 3328 3689 4133 4932 5097 5397 6205 6374 6697 8024

• Agora que temos uma série temporal, vamos plotar:

forecast::autoplot(ts_vendas)

• Método da média simples: Leva em consideração todos os períodos;
• Método “Naive”: Leva em consideração apenas o último período;
• Método da Suavização Exponencial: Leva em consideração a demanda observada nos períodos anteriores, porém faz com que o impacto dos períodos anteriores caia progressivamente. Notação da Previsão: \(\hat{y}_{t+h|t}\)

Previsão para a demanda \(y_{t+h}\), considerando dados disponiveis até \(y_{t}\)

Equação da Previsão: \[\hat{y}_{t+h|t} = \alpha y_{t} + \alpha(1-\alpha)y_{t-1} + \alpha(1-\alpha)^2y_{t-2} + ... \ para \ 0 \leq \alpha \leq 1\]

• Exemplo da Aplicação da Equação:

\[\hat{y}_{t+h|t} = \alpha y_{t} + \alpha(1-\alpha)y_{t-1} + \alpha(1-\alpha)^2y_{t-2} + ... \ para \ 0 \leq \alpha \leq 1\]

Previsão: \[\hat{y}_{t+h|t} = l_t\] Nível: \[l_t = \alpha y_t + (1-\alpha)l_{t-1}\] O valor \(\alpha\) é obtido minimizando os erros ao rodar o modelo em um set de teste.

• A função ‘ses’ (Simple Exponetial Smoothing) estima um modelo de suavização exponencial simples:

dados_petroleo = window(oil, start= 1996)
previsao_petroleo = ses(dados_petroleo, h = 5) # Previsao para os próximos 5 anos
summary(previsao_petroleo)

## 
## Forecast method: Simple exponential smoothing
## 
## Model Information:
## Simple exponential smoothing 
## 
## Call:
##  ses(y = dados_petroleo, h = 5) 
## 
##   Smoothing parameters:
##     alpha = 0.8339 
## 
##   Initial states:
##     l = 446.5868 
## 
##   sigma:  29.8282
## 
##      AIC     AICc      BIC 
## 178.1430 179.8573 180.8141 
## 
## Error measures:
##                    ME     RMSE     MAE      MPE     MAPE      MASE
## Training set 6.401975 28.12234 22.2587 1.097574 4.610635 0.9256774
##                     ACF1
## Training set -0.03377748
## 
## Forecasts:
##      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2014       542.6806 504.4541 580.9070 484.2183 601.1429
## 2015       542.6806 492.9073 592.4539 466.5589 618.8023
## 2016       542.6806 483.5747 601.7864 452.2860 633.0752
## 2017       542.6806 475.5269 609.8343 439.9778 645.3834
## 2018       542.6806 468.3452 617.0159 428.9945 656.3667

• Visualizando a Previsão:

autoplot(previsao_petroleo)

Previsão: \[\hat{y}_{t+h|t} = l_t + h b_t\] Nível: \[l_t = \alpha y_t + (1-\alpha)(l_{t-1}+b_{t-1})\] Tendência: \[b_t = \beta *(l_t-l_{t-1}) + (1-\beta)b_{t-1}\]

• A função ‘holt’ (Simple Exponetial Smoothing) estima o modelo de Holt (que considera a tendência):

previsao_petroleo_holt = holt(dados_petroleo, h = 5) # Previsao para os próximos 5 anos
autoplot(previsao_petroleo_holt)

• Também podemos observar os parâmetros estimados para o modelo:

summary(previsao_petroleo_holt)

## 
## Forecast method: Holt's method
## 
## Model Information:
## Holt's method 
## 
## Call:
##  holt(y = dados_petroleo, h = 5) 
## 
##   Smoothing parameters:
##     alpha = 1e-04 
##     beta  = 1e-04 
## 
##   Initial states:
##     l = 428.4833 
##     b = 5.5771 
## 
##   sigma:  27.8684
## 
##      AIC     AICc      BIC 
## 177.2927 182.2927 181.7446 
## 
## Error measures:
##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -0.2851768 24.57757 20.53231 -0.3285466 4.337459 0.8538816
##                   ACF1
## Training set 0.3429178
## 
## Forecasts:
##      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2014       534.4382 498.7234 570.1529 479.8172 589.0591
## 2015       540.0148 504.3000 575.7295 485.3938 594.6357
## 2016       545.5913 509.8766 581.3061 490.9704 600.2123
## 2017       551.1679 515.4532 586.8827 496.5469 605.7889
## 2018       556.7445 521.0298 592.4592 502.1235 611.3655

• Podemos utilizar o modelo Holt com “Damp”, pressupondo que o crescimento não será linear no longo prazo.

Previsão: \[\hat{y}_{t+h|t} = l_t + (\phi+\phi^2+...+\phi^h) b_t\]

Nível: \[l_t = \alpha y_t + (1-\alpha)(l_{t-1}+\phi b_{t-1})\]

Tendência: \[b_t = \beta *(l_t-l_{t-1}) + (1-\beta) \phi b_{t-1}\] O parâmetro \(\phi\) é entre 0 e 1. Se o parâmetro for igual a 1, o crescimento será linear.

previsao_petroleo_holt = holt(dados_petroleo, h = 15, PI = F)
previsao_petroleo_holt_damp = holt(dados_petroleo, h = 15, damped = T, PI = F) # Previsao para os próximos 5 anos
autoplot(dados_petroleo) + 
  autolayer(previsao_petroleo_holt, series="Holt Linear") + 
  autolayer(previsao_petroleo_holt_damp, series="Holt com Damp")

• Método Aditivo: Adequado quando a amplitude dos ciclos de sazonalidade não está correlacionada ao tempo.
• Método Multiplicativo: Adequado quando a amplitude dos ciclos de sazonalidade está correlacionadao ao tempo.

Previsão: \[\hat{y}_{t+h|t} = l_t + hb_t + s_{t-m+h_m^+}\]

Nível: \[l_t = \alpha (y_t-s_{t-m}) + (1-\alpha)(l_{t-1}+b_{t-1})\]

Tendência: \[b_t = \beta *(l_t-l_{t-1}) + (1-\beta) b_{t-1}\]

Sazonalidade: \[b_t = \gamma(y_t - l_{t-1 - b_{t-1}}) + (1-\gamma)s_{t-m} \] \(s_{t-m+h_m^+}\): componente de sazonalidade do último ano de dados disponíveis. \(m\): Período de sazonalidade. A média do componente de sazonalidade tende a zero.

Previsão: \[\hat{y}_{t+h|t} = (l_t + hb_t) s_{t-m+h_m^+}\]

Nível: \[l_t = \alpha \frac{y_t}{s_{t-m}} + (1-\alpha)(l_{t-1}+b_{t-1})\]

Tendência: \[b_t = \beta *(l_t-l_{t-1}) + (1-\beta) b_{t-1}\]

Sazonalidade: \[b_t = \gamma \frac{y_t}{l_{t-1} + b_{t-1}} + (1-\gamma)s_{t-m} \] \(s_{t-m+h_m^+}\): componente de sazonalidade do último ano de dados disponíveis. \(m\): Período de sazonalidade. A média do componente de sazonalidade tende a 1.

Relembrando a série de anti-glicêmicos.

# Produção de Gás na Austrália
autoplot(a10) + ylab("Demanda AntiDiabeticos")

Relembrando a série de anti-glicêmicos.

# Produção de Gás na Austrália
previsao_anti_diab_aditivo = hw(a10, seasonal = "additive", PI= F)
previsao_anti_diab_multiplicativo = hw(a10, seasonal = "multiplicative", PI = F)
autoplot(a10) + ylab("Demanda de Remédios") +
  autolayer(previsao_anti_diab_aditivo, series="HW Add.") + 
  autolayer(previsao_anti_diab_multiplicativo, series="HW Mult.")

• Colete dados de demanda de diferentes produts em sua empresa;
• Verifique quais são as técnicas de previsão aplicadas a estes produtos;
• Identifique se há ou não tendência e sazonalidade;
• Aplique os métodos aprendidos;
• Proponha um modelo de previsão novo;
• Compare a acurácia do novo modelo em relação ao que existe atualmente na empresa.

• O modelo “naive” simplesmente pressupõe que o futuro repetirá o passado, logo a previsão é correspondente ao último valor observado.

previsao = naive(oil)
autoplot(oil, series="Dados") + 
  xlab("Ano") +
  autolayer(fitted(previsao), series = "Previsao") + 
  ggtitle("Produção de Petróleo na Arábia Saudita")

• Os erros também podem ser plotados. Espera-se que os erros sejam “white noise”. Se isso é verdade, é provável que o modelo tenha captado toda a informação disponível nos dados.

autoplot(residuals(previsao))

• Os erros deveriam ser não-correlacionados;
• Os erros devem ter média zero; Propriedades úteis:
• Variância Constante;
• São normalmente distribuídos.

• A função ‘checkresiduals’ verifica os pressupostos indicados anteriormente. Espera-se que o resultado do p-valor do Ljung-box seja acima de 0.05, indicando que os erros não são correlacionados.

checkresiduals(previsao)

## 
##  Ljung-Box test
## 
## data:  Residuals from Naive method
## Q* = 11.814, df = 9.8, p-value = 0.2824
## 
## Model df: 0.   Total lags used: 9.8

• Ainda assim, as previsões podem ser boas e podem ser utilizadas, porém os intervalos de predição podem ser muito justos ou amplos em função deste problema.

• Traning Set: Parte dos dados que você utiliza para construir o modelo.
• Test Set: Parte dos dados que você utiliza para testar o modelo.
• O test set não pode ser usado para calcular a previsão.
• Um modelo que se ajusta bem aos dados de treinamento não necessáriamente terá uma boa previsão;
• É comum construirmos um modelo altamente complexo que possui poucos erros nos dados de treinamento, porém gera péssimas previsões. Isto é chamado de overfitting, e deve ser evitado.

• Mean Absolute Error: \(MAE = avg(|e_t|)\)
• Mean Square Error: \(MAE = avg(|e_t^2|)\)
• Mean Absolute Percentage Error: \(MAPE = 100 * avg(|e_t/y_t|)\)
• Mean Absolute Scaled Error: \(MASE = MAE / Q\)

• Separamos o modelo em duas partes, e geramos um modelo naive para ilustrar este ponto:

treinamento = window(oil, end=2003)
teste = window(oil, start= 2004)
previsao = naive(treinamento,h = 10) # h = Número de períodos a prever
autoplot(previsao) + 
  ylab("Vendade Petróleo") + 
  autolayer(teste, series = "Dados de Teste")

• Usar o Comando Accuracy para observar o modelo:

accuracy(previsao, teste)

##                    ME     RMSE      MAE      MPE      MAPE      MASE
## Training set  9.87358 52.56156 39.42504 2.506565 12.570647 1.0000000
## Test set     21.60250 35.09832 29.97666 3.963914  5.777875 0.7603458
##                   ACF1 Theil's U
## Training set 0.1801528        NA
## Test set     0.4029519  1.184862

• O agregador de dados Quandl possui milhares de séries temporais disponíveis diretamente no R, pela biblioteca Quandl.
• Exemplo: obtendo a série do índice Bovespa

library(Quandl)
ts_bovespa = Quandl("BCB/7", type = "ts")
autoplot(ts_bovespa)